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We study the flow of energy between a harmonic oscillator (HO) and an external environment consisting 
of N two-degrees of freedom non-linear oscillators, ranging from integrable to chaotic according to a control 
parameter. The coupling between the HO and the environment is bilinear in the coordinates and scales with 
system size as 1 / y/N. We study the conditions for energy dissipation and thermalization as a function of N 
and of the dynamical regime of the non-linear oscillators. The study is classical and based on single realization 
of the dynamics, as opposed to ensemble averages over many realizations. We find that dissipation occurs in 
the chaotic regime for a fairly small N, leading to the thermalization of the HO and environment a Boltzmann 
distribution of energies for a well defined temperature. We develop a simple analytical treatment, based on the 
linear response theory, that justifies the coupling scaling and reproduces the numerical simulations when the 
environment is in the chaotic regime. 



I. INTRODUCTION 

Dissipation is usually modeled by coupling a central sys- 
tem to an external environment. The environment can be rep- 
resented by an infinite collection of modes (either harmonic 
or spin-like) at a given temperature, or by few chaotic degrees 
of freedom. Models in the first category include the Caldeira- 
Leggett (CL) U], where the environment is represented by N 
harmonic oscillators with a linear distribution of frequencies. 
In the thermodynamic limit N — » 00 , the central system obeys 
the classical Langevin equation, exhibiting exponential dissi- 
pation subjected to fluctuating forces. 

In the case of small environments, the exchange of energy 
with the central system exhibits large fluctuations and equili- 
bration is achieved only after averaging over an ensemble of 
realizations, usually of the micro-canonical type. Wilkinson 
10], and later Berry and Robbins [3], studied the coupling of 
a fast chaotic environment to a slow central system using the 
properties of adiabatic invariants. They showed that the adi- 
abatic invariance of the chaotic energy shell [4-6] leads, in 
first order, to a Born-Oppenheimer conservative type of re- 
action on the central system. Second order corrections lead 
to a force proportional to the velocity whose symmetric part 
corresponds to friction and whose antisymmetric part acts as 
a geometric magnetism [0 0] . Since the publication of these 
pioneering works, a large number of papers has explored the 
possibility of dissipation caused by small environments, both 
classically ll7l-[TTIl and quantum mechanically E3-H2], where 
decoherence is the main focus. 

Although the approach to dissipation via ensemble average 
over small environments has lead to several interesting results, 
it cannot account for situations where a single experiment is 
performed and no averages take place. This is, for example, 
the case of the C^q molecule interacting with its own internal 
degrees of freedom. As pointed out in [23], couplings of these 
molecules to the environment or to their own internal degrees 
of freedom play an essential role in the appearance of decoher- 
ence, destroying self-interferences in a double slit experiment, 
and deserve a careful analysis. When the environment has rel- 
atively few degrees of freedom it is not clear whether they 
lead to standard dissipation effects over the collective degrees 
of freedom. 



The main goal of this paper is to study dissipation from 
the classical point of view as induced by small environments 
without resorting to ensemble averages. To distinguish this 
process from those obtained after averaging, we term it ef- 
fective dissipation. We follow the general approach of Feyn- 
man and Vernon l24ll where the system of interest, or central 
system, has much fewer degrees of freedom than the environ- 
ment, which, however, is also small. In our simulations the 
central system is a harmonic oscillator with a single degree of 
freedom whereas the environment has 2N freedoms. We study 
the system as a function of N and of its dynamical character, 
regular or chaotic. We wish to determine the minimum value 
of N for which the environment behaves as a thermal bath, 
inducing irreversible dissipation on the harmonic oscillator. 



II. MODEL 

We consider a globally conservative system governed by 
the Hamiltonian 



H = Hho +He + X^Vi 



where 
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The environment He is modeled by the finite sum of iden- 
tical Quartic Systems (QS) [25]. Each QS has two degrees 
of freedom and is invariant under a scaling transformation in 
such a way that the dynamics in all energy shells are the same. 
This property allows us to adjust the QS's energy and time 
scales without changing its dynamical regime, which is solely 
controlled by the parameter a. The QS is integrable for = 1.0, 
strongly chaotic for a < 0.1 and mixed for intermediate val- 
ues. 

In addition to the dynamical regime, the number of QS's 
coupled to the harmonic oscillator (HO) plays a crucial role 
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in the dynamics of dissipation. In order to compare results 
for different values of N, the coupling parameter Xn, as seen 
in equation (Q]), has to be properly normalized with the size 
of the environment, so that the effective coupling becomes 
independent of N for large N. We use Xn = X / y/N and derive 
this scaling later using linear response theory. 

We want to study the behavior of the system for a single 
global initial condition and we consider two types of initializa- 
tion. In both cases the HO starts from q = and p = y/2mE Q 
with m = 1, (0 o = 0.3 and E = 10.0. In the first type, the 
initial condition for each QS is chosen randomly at a fixed 
energy shell Eq. This choice is related to the microcanonical 
distribution 

1 N 

P=?U 5 ( H Cn-Ecn), (5) 

where Ec n = Ec, for n = 1 , . . . , N. We refer to it as pseudo- 
microcanonicaU to distinguish it from the real microcanonical 
distribution which is uniform over the entire energy shell of 
the Af QS's. 

The second type of initialization corresponds to pick a dif- 
ferent energy shell Ec n for each QS from the Boltzmann dis- 
tribution p(E) = exp (-E I Eqs) I Eqs and we call it pseudo- 
canonical. The initial condition for each QS is again chosen 
randomly at the corresponding energy Ec n - 



III. NUMERICAL RESULTS 

In all simulations the coupling parameter will be fixed at 
X = 0.01. We consider initially the pseudo-microcanonical 
distribution with Ec = 0.01. 

We first study the role of chaos in promoting dissipation. 
Fixing N = 100 we integrate Hamilton's equations for differ- 
ent values of the parameter a and compute the energy of the 
harmonic oscillator Hq as a function of the time. Time is mea- 
sured in units of the HO's period, t = 2n/co « 20.94. FigQ] 
shows the results of simulations for a single numerical real- 
ization for three values of a. 

In the integrable regime, a = 1.0, the effect of the environ- 
ment is minimal, causing a small decrease in the oscillator's 
energy and inducing small amplitude, noisy, oscillations. In 
the mixed case, a = 0.5, the energy loss is slightly larger and 
the oscillations are more clearly visible, particularly at shorter 
times. At long times the energy fluctuates around 0.8E o . We 
refer to the papers [[Tol Qj]] for numerical results on analogous 
systems, computed via ensemble averages, in the regular and 
in mixed regimes. The authors describe in great detail the de- 
pendence of the dissipation rate on the number of degrees of 
freedom of the environment in these regimes. 

In the chaotic regime the energy of the HO shows a clear 
exponential decay for N = 100 QS's. Similar results were ob- 
tained in [9] with N = 1 but performing ensemble averages 
over several thousands of trajectories. The results displayed 
in FigQ] are for a single trajectory. Changing the initial con- 
ditions in the environment, but keeping the same microcanon- 
ical constraint, produces energy curves that differ only in the 




20 40 60 80 100 

t/x 



Figure 1. (Color online) Energy of the HO when coupled to N = 100 
QS in distinct dynamical regimes: integrable, a =1.0 (blue); mixed, 
a = 0.5 (red); chaotic, a = 0.01 (black). 

details of the fluctuations, but not in their main qualitative fea- 
tures. 

The number of degrees of freedom of the environment is of 
fundamental importance to the effective dissipative behavior. 
The exponential decay observed in FigQ] occurs only if N is 
sufficiently large. Figj2] shows the typical behavior of Ho(t) 
for different values of N. 
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Figure 2. Dissipation curve as function of the number of QS in the 
environment. 

For small environments QS (1 < N < 6), the fluctuations 
are so large that dissipation cannot be characterized for each 
numerical realization. In these cases, and in particular for 
N =1, exponential dissipation can be observed only in an en- 
semble averaged treatment of the energy H. 

Increasing the number of QS's to 7 < N < 20 makes the 
decreasing behavior of the HO energy more evident, although 
the fluctuations are still large. Finally, for the parameters of 
this simulation, the curves converge to a ^-independent ex- 
ponential decay H Q (t) ~ e' 71 for TV > Afc — 100. Such con- 
vergence occurs only if the coupling constant scales properly 



3 



with N. Clearly Xn has to decrease with the number of QS 
in the environment, otherwise the interaction Hamiltonian Hi 
would diverge as N — » ©o. Only if Xn = X/ vW, the effect of 
the environment over the HO converges as N becomes large, 
as shown in Figj2j This scaling of the coupling constant will 
be demonstrated using linear response theory. 

Environments with N > Nq behave as infinite reservoirs. 
Clearly, Nq depends on the model parameters and, in particu- 
lar, on a, that controls the degree of chaos in the QS. Indeed, 
dissipation is intimately associated with correlation functions 
that contain information about the memory of the environ- 
ment. Increasing the phase space region covered by chaos (by 
decreasing a) makes these correlation functions decay faster, 
leading to a faster transfer of energy from the oscillator to 
the environment, even for small values of N. This makes Nq 
smaller for smaller values of a. 

The flow of energy between the HO and environment also 
depends on the relative time scales of the HO and QS. Increas- 
ing the HO frequency (co ) decreases the relative velocity of 
the two subsystems, leading to a less efficient energy dissipa- 
tion. However, since the QS is scalable, we can always change 
its energy in such a way to readjust its relative velocity. 

The decay rate 7, or the characteristic time t^s = 1/7' ^ so 
changes with A, as shown in Fig J5] for N = 100 and a = 0.01. 
The dependence is well fitted by a quadratic law, a result that 
will also be demonstrated using linear response theory. 
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Figure 3. Inverse characteristic time (7) as a function of the coupling 
constant (A) for N = 100 and a = 0.01. The curve is well fitted by a 
quadratic function. 

The numerical results in Fig J2] show that small chaotic envi- 
ronments may behave as infinite reservoirs, absorbing energy 
from the HO. This behavior can be further characterized by 
associating a temperature to the full system after its equilibra- 
tion. In order to do so, we look at the energy distribution of 
the HO and the environment at equilibrium. This can be done 
in two different ways: the first is to count, at a fixed time, the 
number of QS's with energy in the interval (E,E + dE) for 
a large number of initial conditions (ensemble average). The 
second possibility is to count, for a single realization, the num- 
ber os QS's with energy in the interval (E,E + dE), starting 



at some time after equilibration and repeating the procedure 
several times at a pre-specified time step (time average). The 
same procedures can be applied to the energy of the HO. The 
statistical distribution of these results will be the same if the 
system is ergodic. We performed numerical calculations using 
both the time and ensemble averages and found very similar 
results, indicating that both the environment and the HO are 
indeed ergodic in equilibrium. 

In order to define a temperature for the environment we use 
the general equipartition theorem [12611 to show that Eqs = 
\Ero- The energy distribution of the environment for N = 
100 is well fitted by exp(-E/E QS ) with E QS « 0.1094. The 
value of Eqs is a measure of the environment's temperature. 
Fig shows that the energy distribution of the HO also follows 
a Boltzmann-like distribution, 

Poh(E) = -J- e - E / £ *o (6) 

where E H q = 0.0723 w 2/3Eq$, showing that the systems are 
indeed in thermal equilibrium. 
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Figure 4. Energy distribution for HO with N = 100 QS's. The tem- 
perature is E H0 = 0.0723 ± 0.0002. Energy was sampled 100,000 
times at equal intervals. 

This is an important result of this paper, that distinguishes 
between the coupling of the HO with a single QS (where re- 
sults are averaged over many realizations with different initial 
conditions) [ 9] and the coupling with many (N ~ 100) simul- 
taneous QS (and no average over realizations). In the former 
case a temperature can be defined only by modifying the defi- 
nition of entropy |27l|28|]. When several QS's are present, the 
indirect, and therefore weak, interaction between the members 
of the environment plays a key role in the dynamics of the sys- 
tem. It is this weak interaction that allows for their eventual 
equilibration. 

These numerical results suggest that classically irreversible 
processes may indeed occur in nature if the system interacts 
with a small, but not too small, chaotic environment. The pri- 
mary coupling between the system and the environment leads 
to dissipation whereas the secondary interactions between the 
few, but chaotic, degrees of freedom of the environment, leads 
to thermalization. 



4 



IV. LINEAR RESPONSE 

We can understand the basic mechanisms of dissipation us- 
ing linear response theory (LRT) 1 2911 . We start by writing an 
approximation to the full Hamiltonian of the system as 

H = H E (Q,P)+Hj{Q,P,t) 

where Q and P represent the full set of canonical variables of 
the environment, {x n ^y n ^Px n ^Pyn)-> n ~ 1 , • • • . The second 
term Hi(Q,P,t) = A(Q,P)x(t) is a perturbation that includes 
the central system through the time dependent function x(t)- 
This approximation to (03 is valid in the limit of weak cou- 
pling, since the central system is treated as an external source 
that does not respond to the environment. The feedback of 
the environment on the central system will be considered in a 
moment. 

The dynamics can be described by the Liouville equation 



^ =iL E p(t) + iL I (t)p(t). 

where the Liouville operators L E and L/ are given by 
dH dp dH dp 



(7) 



iLp = £ 



dQj dPj dPj dQj 



{H,p}. (8) 



with H replaced by H E and Hi respectively. Eq.© can also 
be written in integral form as 

p(t) = e^-^p^ + i f e^- s ^Li(s)p(s)ds. (9) 

JtQ 

The initial distribution p(*o) is assumed to be invariant 
under H E , like, for example, the pseudo-microcanonical or 
pseudo-canonical distributions given by Eq.©. The action of 
Hi(t) removes the environment from this initial equilibrium. 
If the perturbation is small, Eq.© can be expanded to first 
order in L/ as 



p(t)=p(t )- 



JtQ 



s)Lq 



{Hi(s),p(t )}ds. (10) 



The ensemble average of a general function B(Q,P) can then 
be calculated as 

(B(Q,P))(t) = (B(Q(t),P(t)))o+ f foA(t-s)x(s)dQdPds. 

JtQ 

(ID 

where 



$BA(t) = ({A(Q,P),B(Q(t),P(t))}) c 



(12) 



is the response function of B(Q, P) when the environment is 
under the influence of the perturbation A (g, P). The subscript 
indicates that the averages are computed with the initial in- 
variant distribution p (Q, P, to) • 

We can now compute the response of the central system to 
the environment. For the Hamiltonian (OQ) the HO satisfies the 
equation 



and is perturbed by the 'external force' F(t) = ^-X(t). If the 
dynamics of the bath coordinates x n (t) are chaotic we may 
replace X(t) by its average (X(t)): 



m 



(14) 



Using Eq.lfTTT). 

(X(t)) = (Yxn(t)) -hi f<Pxx(t-s)q(s)ds. (15) 

where ®xx(t-s) = {{X(t),X(s)}} . The first term in the 
righthand side of Eq.(fT3T> is zero due to the parity of He- Using 
the pseudo-microcanonical as the initial invariant distribution 
it follows that 



&xx(t-s)= Y t <t£\t-s)=N* 



(16) 



n=\ 



where 



®xx(t-s) = 



4E C (0) 
t-s 



4£ c (0) 



(x(t)P x (s)) + 
(P x (t)P x (s)) , 



(17) 



is the Response Function obtained by Bonanea |9] for the case 
of a single QS coupled with the oscillator. Substituting dT6l> in 
(ITSb and ([HI) we obtain 



(x>lq : 



X 2 N N 



f 
Jo 



<S> xx (t-s)q(s)ds, 



(18) 



that shows that it is necessary to re-scale the coupling constant 
as An = h/y/N, in order for the external driving F(t) to be 
independent of the size of the environment in the limit of large 
N. 

The approximation where the N chaotic systems forming 
the environment are treated independently correspond to en- 
semble averages of a single QS coupled to the central system. 
As we have seen, this approximation cannot describe the long 
time equilibration of the environment, since this requires some 
interaction between the chaotic systems. Linear response the- 
ory cannot, therefore describe the long-time behavior of the 
system. In particular, it cannot describe the change from an 
initial pseudo-microcanonical distribution to the final canoni- 
cal distribution. 

Fig(5] shows the breakdown of the pseudo-microcanonical 
distribution and its consecutive rearrangement to a 
Boltzmann-like distribution as a function of time. The 
transition between the distributions occurs at about 1/10 
of the characteristic time for dissipation, as can be seen 
from Figs 13 and [TO] As pointed out in the last section, 
this redistribution of energy is promoted by the secondary 
interactions between the QS's. 

This rapid equilibration on the environment hints to the pos- 
sibility of using LRT to describe the long time behavior of the 
system by changing from the initial pseudo-microcanonical 
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Figure 5. Time evolution of the environment energy distribution. At t=0 p# is given by equation with Ec n = 0.01 for all QS's and N = 100. 
Note the different scales in the panels. 
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Figure 6. Energy distribution of the environment at different times 
for N = 100. At t = p(E) = exp (-E/E QS ) /E QS with E QS = 0.01. 



to the pseudo-canonical distribution. If the initial tempera- 
ture, as given by /5qs = Eqs, does not change much during the 
time evolution, the environment will not be largely affected 
by the central system, and the interaction between the QS's 
can be ignored. Fig J5] shows an example with N = 100 and 
Eqs (0) = 0.01 where the distribution of energy of the envi- 
ronment does not change significantly with time. 

Therefore, for N sufficiently large, we can establish a 
regime where the central system dissipates energy exponen- 
tially, but the environment's temperature remains nearly con- 
stant. This allows us to use the LRT to describe the system 
dynamics for much longer times. In order to do this we need 



to obtain explicit expressions for Eqs.(fT6l and (fT71> for the 
pseudo-canonical distribution. We write the correlation func- 
tions in dT71) as 



and 



®xx(t- 



(19) 
(20) 



(t-s) d 2 



4 dtds \r^(o) 



-x n (t)x n {s) ) (21) 



The individual autocorrelation functions c n (t,s) 
(x n (t)x n (s)) were obtained numerically and follow the 



typical behavior expected for chaotic systems [30] which can 



be fitted by c n (t,s) = o n e 



-Ut-s) 



cos[v n (t — s)}. However, the 



total correlation function, C(t,s) = \Ym e (o) x n(?) x n{s) 
smooths out the oscillatory behavior exhibited by a single 
QS, as shown in FigJTJ It is important to note that the time 
scale where C(t , s) is significant is much smaller than the time 
scale of the dissipation, as can be seen by comparison with 
FigE 

Therefore, in the time scale of dissipation, we can make the 
approximation 
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Figure 7. The total correlation function C(t, 0) = c w (f , 0) with w = 
l,...,Af = 100 QS's obeying a pseudo-canonical distribution with 
% = 0.01 



where the parameter k e will be adjusted later and = jU^N 
is the maximal amplitude of Cjy(M)- It turns out that A^y in- 
creases linearly with N, as shown in Fig 03 and we call fl E the 
slope of such linear function. The index "E" on K E and \i E 
expresses the dependence with the mean energy of reservoir, 
as we will see in a moment. 

Equation (l22b indicates that the presence of chaos intro- 
duces a fast memory loss in the microscopic dynamics of 
reservoir, so that its dynamics can be described by a Marko- 
vian process, unlike the integrable and mixed regimes. 



Substituting (f)xx(t — s) in equations (ITU) and (fT5l) , and com- 
puting the integrals, we obtain 



CQ 2 q- 



-7M = 



where 



Je 



lX 2 K E ji E N 



lX 2 K E \i E 



(24) 



(25) 



4m 4m 
The inverse characteristic time y E obtained from LRT de- 
creases exponentially with E as shown in Fig |9] (circles). The 
numerically obtained values, on the other hand, displays a 
more complex behavior, as shown by the squares in the same 
figure. The numerical results can be divided into two regimes: 
the first, where < E < 0.01 and y(E) increases with E and, 
the second regime, where E > 0.01 and y(E) decreases with 
E. 
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Figure 8. Linear dependence between the maximal amplitude of the 
total correlation function C(t , 0) and the number of QS's in the envi- 
ronment. The slope depends on the environment's mean energy. 

Using (f22l) we obtain the response function 
5k e \i e N d 



K E n E N(t-s) d 2 
4E Cn (0) dtds 



8(t-s). (23) 



Figure 9. The inverse characteristic time as a function of the envi- 
ronment's mean energy for: theoretical model (y E , red circles); and 
numerical data (7, black squares). 

The discrepancies between the numerical and LRT results 
can be understood if we note that in the first regime, where 
E is small, the QS's are very slow. Since the value of Ye i s 
related to the maximal amplitude of the correlation C(t,s), 
this value increases considerably in these regimes because 
changes in the state of the quartic oscillators are slow. This, 
in turn, implies a super-estimated value for y E . In this regime, 
as E increases, the influence of the environment on the HO 
increases, facilitating the flow of energy and increasing y E . 
This is a non-trivial effect, since an environment with lower 
mean energy should correspond to a cooler thermal bath, with 
a greater dissipation rate. This shows that the finite set of QS's 
cannot be traded by a thermal bath at very low energies. 

After a threshold is passed /recovers the expected decreas- 
ing behavior and the association between the mean energy and 
temperature is restored. Even in this case, where the qualita- 
tive behavior of y E and 7 is the same, the numerical value of 
the coefficients needs some correction. We define the param- 
eter k e = 7/ y E in order to adjust the theoretical model with 
numerical results. When this is done the numerical and theo- 
retical curves fall on top of each other, as shown in FigJTOl 
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Figure 10. Comparison between the numerical results (continuum 
black line) and theoretical model (dashed red lines), for N = 100 and 
k e = 3.68. 



V. CONCLUSIONS 

We constructed a simple model of classical dissipation 
where a harmonic oscillator interacts with an environment 
consisting of N identical quartic sub-systems (QS) with two 
degrees of freedom. Interactions among the sub-systems oc- 
cur via the harmonic oscillator and are, therefore, of second 
order in the coupling constant. In order to compare results 
with different number of QS's we re-scaled the coupling con- 



stant according with Xn = A / y/N. 

The main result of this paper is that a small chaotic envi- 
ronment can indeed behave as an effective infinite reservoir, 
promoting energy flow and equilibration with smaller systems 
in a single realization of the dynamics. This is to be con- 
trasted with previous results where a single chaotic system 
plays the role of environment and thermodynamical behavior 
is achieved via ensemble average over many realizations with 
random initial conditions [9] . For the same values of parame- 
ters the results with N = 100 are roughly equivalent to averag- 
ing over 30,000 initial conditions. Comparing our model en- 
vironment with those composed of harmonic oscillators with 
linear distribution of frequencies [OH [Till , the number of os- 
cillators needed to mimic an infinite reservoir is at least an 
order of magnitude larger than needed for chaotic QS's. 

At long times the environment and the HO equilibrate and 
their energy distributions converge to Boltzmann-like curves 
with the same temperature. The equilibration of the environ- 
ment depends on the indirect interactions among its members 
that occur via the HO. 

We have also developed a Linear Response Theory to 
describe the system. We have shown, in particular, that the 
scaling of the coupling constant with N is indeed correct and 
that dissipation constant y depends quadratically on A . 
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